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^ ■ Abstract 

Q . A simple reweighting scheme is proposed for Monte Carlo simulations of in- 

teracting particle systems, permitting one to study various parameter values 
in a single study, and improving efficiency by an order of magnitude. Unlike 
earlier reweighting schemes, the present approach does not require knowl- 
edge of the stationary probability distribution, and so is applicable out of 
C^ ■ equilibrium. The method is applied to the contact process in two and three 

C . dimensions, yielding the critical parameter and spreading exponents to un- 

rj^ \ precedented precision. 
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Reweighting and histogram methods, which permit one to calculate thermal averages at 
different temperatures from a single sample, have greatly improved the efficiency of Monte 
Carlo simulations [0]. Given a sample of configurations, Ci,...,CAr, generated at temper- 
ature T, one can in principle generate a sample appropriate to some temperature T', by 
assigning the weight Wi = Pi{T') / Pi{T) to configuration i. Until now, such methods have 
been restricted to equilibrium, where the stationary probability distribution is Pi oc e~^'/'^^. 
Away from equilibrium, one does not in general know the stationary distribution, so there 
is no simple way to evaluate the Wi. Recently Grassberger and Zhang showed how a "self- 
organized" formulation of directed percolation can be used to study several parameter values 
in a single run, without reweighting 0. In this work I introduce a reweighting scheme for 
interacting particle systems, based on the observation that, despite our ignorance of the 
stationary distribution on configuration space, we can write down the probability for any 
sequence of events between time zero and time t. 

I apply the reweighting method to the contact process (CP), a simple particle system 
(lattice Markov process) exhibiting a phase transition to an absorbing (frozen) state at 
a critical value of the creation rate |^. This model belongs to the universality class of 
directed percolation [Q and Reggeon field theory p, (it is one of the most well-studied 
representatives of this class [Q), and is pertinent to models of epidemics |^, catalysis [Q, 
and damage spreading |^ . In the CP each site of the hypercubic lattice Z"^ is either vacant 
or occupied by a particle. Particles are created at vacant sites at rate Xn/2d, where n is the 
number of occupied nearest-neighbors, and are annihilated at unit rate, independent of the 
surrounding configuration. The order parameter is the particle density p; the vacuum state, 
p = is absorbing. As A is increased beyond Ac, there is a continuous phase transition from 
the vacuum to an active state; for A > Ac, p ~ (A — Ac)^ in the stationary state. 

There are a number of ways (equivalent as regards scaling behavior), of implementing the 
CP in a simulation algorithm; in this work I follow the widely used practice of maintaining 
a list of all occupied sites. Trials begin at time zero, from a fixed initial configuration. 
Subsequent events involve selecting (at random) an occupied site x from the Np sites on 
the list, selecting a process — creation with probability p = A/(l + A), annihilation with 
probability 1 — p — and, in the case of creation, selecting one of the 2d nearest neighbors, 
y, of X. (The creation attempt succeeds if y is vacant). The time increment At associated 
with an event is 1/Np, where Np is the number of occupied sites immediately prior to the 
event. A trial ends when all the particles have vanished, or at the first event with time > tm, 
a predetermined maximum time. 

Consider a single trial, extending from time zero up to t^; for simplicity, suppose that 
initially there is but a single particle, located at the origin. (If all of the particles disappear 
at some time t', the trial is trapped in the vacuum state for all later times.) A trial consists 
of sequence S of events, each involving the annilihation or (attempted) creation of a particle. 
With the help of diagrams reminiscent of the "percolation substructure" invoked in defining 
the CP ||TO|JTT[] , we can write down the probability of sequence S; examples are shown in 
Fig. 1. Each annihilation event carries a factor of (1 —p)/Np, each creation event (succesful 
or not) a factor p/{2dNp). The probability Pp{S) of a sequence S, extending to time t, is 
simply the product of all factors associated with events occuring at times t' < t. Now, for 
finite t, the set St of possible sequences is finite, and if A{S; t) is any property of the system 
(e.g., the number of particles at time t), then its expectation is 



(A),= Y. PpiS)A{S;t) . (1) 

S£St 

In a Monte Carlo simulation we generate a sample Si,...,Sn, drawn from the distribution 
Pp{S), which yields the estimate 

1 N 

A;p^j;^T.MSk;t). (2) 

k=l 

From our analysis of Pp{S), it is evident that its p-dependence only involves the numbers 
c{S) and a{S) of creation and annihilation events, respectively; the ratio of the probabilities 
associated with two different values of p is 

Thus the reweighted estimate, 

1 N 

A.y^-Y.^iSk)A{Sk;t), (4) 



has expectation 
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Y.PpiSMS)A{S;t) = {At)p>, (5) 
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and is an unbiased estimator of A in the process with creation probability p'. 

As in other applications of reweighting, it is not enough to have an unbiased estimator; 
one must also ensure that the sample generated with creation probability p has a reasonable 
degree of overlap with a typical sample at p'. We expect that as the duration tm increases, 
the range of A values for which a sample is useful will narrow. To estimate this range of 
values, consider a sequence of r events. The probability that exactly c of these are creation 
attempts is given by the binomial distribution, 

^(^'^) = c!(r-c)!(l + Ar ^^^ 

so that (c) = rp and the standard deviation cr is given by cr^ = rA/(l + A)^. In a typical 
sequence generated with creation probability p, the actual number of creation events will be 
in the range rp ± a, corresponding to a creation rate of A' ~ A[l ± (1 + A)(Ar)^^/^]. Thus 
the typical fluctuation in A is 




5A=(1 + A)W-. (7) 



In the two-dimensional CP, for example, one finds that r ~ 0.7 t^'^^ for A ~ Ac — 1.6488. 
This corresponds to 6X = 0.013 for tm = 1000. Although the range may appear narrow, it 
is more than sufficient for "time-dependent" simulations, which typically focus on a small 
interval near A^.. 



I have applied reweighting to the CP in two and three dimensions, in studies using two 
kinds of initial states. In the case described above, one studies the survival and spread of 
activity from a single 'seed' particle. (The system size must be sufficient to guarantee that 
particles never reach the boundaries.) The principal quantities of interest are the survival 
probability P{t), the mean number of particles n{t), and the mean-square distance R^ of 
particles from the origin. At the critical point these are known to follow asymptotic power- 
laws [|g, 

P{t) ~ r' , (8) 

n{t) ~ t" , (9) 

and 

R'it) ~ f . (10) 

Away from Ac these quantities show deviations from power laws. Since the CP exhibits 
corrections to scaling of the form P(t) ~ t^^[l + at'"' + ■ ■ ■] (similarly for n and R"^), it 
is useful to plot derivatives (local slopes of the log-log plots) versus t~^ when extracting 
the critical exponents 6, rj and z |jl2|. These exponents are connected by the hyperscaling 
relation: 46 + 27] = dz. 

The second kind of study begins with all sites occupied, and follows the decay of the 
particle density p{t). In this case the signature of the critical point is power-law decay, 
p ~ t~^ in the short-time regime, i.e., before the correlation length has attained the system 
size L. A scaling argument implies 6 = 6 |r5| , P |. 



As a preliminary test, I compared two results for the survival prbability P{t) in the 
two-dimensional CP at A = 1.665, one obtained by reweighting a sample of A^ = 10^ trials 
generated at A = 1.650, the other from a similar run but without reweighting (A = 1.665). 
The relative difference between the two results for P{t) remains < 0.017 for t < tm = 
400. Since each result has a relative uncertainty of y (1 — P{t))/NP{t) ~ 0.01 for t = t^, 
{P{tm) — 0.085), the difference between the two results is fully consistent with sample-to- 
sample fluctuations. 

In two dimensions I performed spreading simulations extending to t^ = 2980, on lattices 
of up to 1200 X 1200 sites. Samples generated at a central value, Aq, were reweighted so as to 
study ten additional values, A = Aq ± mAA, with m = 1, ..., 5. The general strategy is to use 
relatively small samples and run times initially, to bracket Ac, and then extend the sample 
size and run time to make finer distinctions. Thus a sample of 10^ trials with tm = 665 and 
AA = 10~^ is already sufficient to restrict Ac to the interval [1.64875, 1.64895]. The most 
sensitive indicator of criticality is the local-slope plot of rjt, defined as the derivative of Inn 
with respect to Int. [Numerically rit is estimated from a least-squares linear fit to the Inn 
data for a set of n^ = 17-25 equally-spaced values (an increment of 0.1) of Int; it is plotted 
versus t~^, ta being the geometric mean of the t- values over the rij intervals. 6t and Zt are 
obtained similarly.] 

To refine the estimate for Ac I generated two samples of 2.5 x 10^ trials each, with 
tm = 2980, Aq = 1.64880, and AA = 4 x 10"^. The critical point was determined from the 
plot of rjt, which shows a clear deviation from smooth behavior for off-critical values. (Fig. 



2 shows the rjt plot for one of the two runs. Note that the leading correction to scaling 
does not follow a simple 1/t decay, as was also noted for DP in 2+1 dimensions [0.) These 
studies yield Ac = 1.64877(3), the number in parentheses denoting the uncertainty in the 
last figure. Extrapolating the local-slope plots, I obtain 6 = 0.4523(10), r] = 0.2293(4), and 
z = 1.1316(4), which are in good agreement with, and generally sharper than, the results of 
previous large-scale simulations of directed percolation 0, and of the ZGB surface catalysis 
model |TB[ (see Table I). 

Grassberger and Zhang p| noted that the derivatives of lnn{t) and InP(t) with respect 
to A (evaluated at Ac), grow ~ t^^'^ii. Analyzing n(t) in this fashion yields i/j| = 1.292(4), in 
good agreement with their estimate of 1.295(6). 

In three dimensions, I performed four runs of 10^ trials each, extending to tm = 2208. (To 
avoid finite-size effects, I do not use an occupancy array in the three dimensional simulations, 
but simply search the particle list to detect overlaps.) Two studies used Aq = 1.31686, 
AA = 3 X 10"^ in the others, Aq = 1.31689 and AA = 2 x 10"^ Fig. 3 shows the local 
slopes for one of the latter runs. While the extrapolation of t] is straightforward, the strong 
linear correction to 6 {6t — 6 — 4.96t^^ + ■ ■ O? renders it advantageous to subtract the linear 
term when estimating the exponent (inset of Fig. 3b). In the case of z, the curves for all A 
values are virtually identical. It is difficult to extrapolate the 1/t plot; the present estimate 
is derived by extrapolating the local slope plotted versus t~^/^, as in Fig. 3c. Averaging over 
the results for the four runs yields the following estimates for the three dimensional CP: 

Ac = 1.31686(1), r/ = 0.110(1), 

S = 0.7263(11), z = 1.042(2) 

(the uncertainties represent one standard deviation). These are in good accord with hyper- 
scaling: 46 + 27] — 3z = —0.001(12). The present results are compatible with Jensen's esti- 
mates of Ac = 1.3168(1), T] = 0.114(4), and 5 = 0.730(4), but not with his value z = 1.052(3) 
p!6| . Analysis of d\nn(t)/d\ yields z/|| = 1.114(4), just consistent with Jensen's result, 
1.105(5). 

I studied the initial decay of the density (starting with all sites occupied), in the two- 
dimensional CP for system sizes L = 32, 64, and 128. Fig. 4 shows the results of one of 
three sets with L = 128, 10^ trials, Aq = 1.6490, and AA = 10"^ (The studies for L = 32 
and 64 used similar parameters, in runs of 5 x 10^ and 2 x 10^ trials, respectively.) For 
the relatively short runs employed here {tm < 1096), all trials survive up to t^- While the 
p{t) curves (inset) for different A values are indistinguishable, the local slopes (main graph) 
vary quite systematically with A. Since the power law p ~ t~^ obtains for finite times not 
t — i> oo, I plot the local slope versus Int in this case. One can distinguish three regimes: an 
initial phase in which St increases for all A values, a late stage in which it decreases, as p{t) 
approaches a value ~ L~^^'^^ as predicted by finite size scaling, and an intermediate regime 
in which 5t is more or less constant. Associating Ac with the 5t most nearly constant in the 
intermediate regime yields Ac = 1.6492(1) for L = 32, 1.6491(1) for L = 64, and 1.64898(5) 
for L = 128. The corresponding estimates for 6 (from the fiat portion of each curve) are 
0.4508(10) for L = 32 and 0.4520(5) for L = 64 and 128. Thus the Ac estimates appear to 
be approaching the value derived from spreading simulations; the two kinds of studies yield 
consistent results for 5. 



In summary, I propose a simple reweighting scheme for nonequilibrium lattice models, 
and apply it to the contact process in two and three dimensions. Since spreading simula- 
tions are usually repeated for five or so different A values, reweighting yields roughly an 
order-of-magnitude speedup. In addition to improving efficiency, using the same sample to 
study all parameter values eliminates the effects of independent fluctuations, which compli- 
cate determination of Ac and the critical exponents. These computational advantages have 
permitted determination of the critical parameters of the three-dimensional CP to unprece- 
dented precision. One may expect reweighting to find wide application in simulations of 
nonequilibrium critical phenomena. 
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TABLES 



TABLE I. Critical exponents for DP in two dimensions. 



exponent 


Ref. II 


Ref. [H 


Pres. work 


5 

7] 

z 

A5 + 2ri- 2z 


0.451(3) 
0.229(3) 
1.133(2) 
-0.004(22) 


0.4505(10) 
0.2295(10) 
1.1325(10) 
-0.004(8) 


0.4523(10) 
0.2293(4) 
1.1316(4) 
0.005(6) 



FIGURE CAPTIONS 

FIG. 1. Examples of event sequences in the one-dimensional contact process starting from 
a single particle. Vertical lines represent particles, whose birth is marked by a dot, and 
annihilation by x. Solid and dashed horizontal lines represent, respectively, successful and 
unsuccessful creation events. Probabilites are listed beneath each sequence. 

FIG. 2. Local slope plot for exponent rj in the two dimensional CP. The middle curve (with 
data points) marks Aq = 1.64880. Curves above and below are for A values at intervals of 
AA = 4 X 10~^ above and below Aq. 

FIG. 3a. Local slope plot for exponent 6 in the three dimensional CP. Inset: detail of 
6 — 4.96t~^. The middle curve (with data points) marks Aq = 1.31689; curves above and 
below are for A values at intervals of AA = 2 x 10~^. 

FIG. 3b. Local slope plot for exponent rj in the three dimensional CP. Symbols as in Fig. 
3a. 



FIG. 3c. Local slope plot for exponent z in the three dimensional CP. The inset shows the 
same data plotted versus t~^/^. 

FIG. 4. Local slope plot for exponent 6 in the two dimensional CP starting with all sites 
occupied, L = 128. The data points mark Aq = 1.6490, AA = 5 x 10^^. Inset: decay of the 
density p. 
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